{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Metrics in GPyTorch\n",
    "\n",
    "In this notebook, we will see how to evaluate GPyTorch models with probabilistic metrics. \n",
    "\n",
    "**Note:** It is encouraged to check the Simple GP Regression notebook first if not done already. We'll reuse most of the code from there.\n",
    "\n",
    "We'll be modeling the function\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "  y &= \\sin(2\\pi x) + \\epsilon \\\\\n",
    "  \\epsilon &\\sim \\mathcal{N}(0, 0.04) \n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we set up the train and test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training data is 100 points in [0,1] inclusive regularly spaced\n",
    "train_x = torch.linspace(0, 1, 100)\n",
    "# True function is sin(2*pi*x) with Gaussian noise\n",
    "train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * math.sqrt(0.04)\n",
    "\n",
    "test_x = torch.linspace(0, 1, 51)\n",
    "test_y = torch.sin(test_x * (2 * math.pi)) + torch.randn(test_x.size()) * math.sqrt(0.04)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we define a simple GP regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will use the simplest form of GP model, exact inference\n",
    "class ExactGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
    "    \n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "# initialize likelihood and model\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = ExactGPModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model is ready for hyperparameter learning, but, first let us check how it performs on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAADGCAYAAADSWvcAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4G0lEQVR4nO3deXiU5bn48e8z+5Z9gRCSsBPCJotYlxZxKypgrdoKtv7cigeOyznIsYu26LG29lRrtdbW5WjrqUXqhorWClWKCqgg+yoBwk72dfaZ5/fHJCEJWckkmST357pyJZl58849Ibz3+2z3o7TWCCGEEIaeDkAIIURskIQghBACkIQghBCiliQEIYQQgCQEIYQQtSQhCCGEAKKQEJRSNqXU50qpLUqpHUqpB6MRmBBCiO6lOrsOQSmlAKfWulopZQY+Ae7WWq+PRoBCCCG6h6mzJ9CRjFJd+6259kNWuwkhRC8TlTEEpZRRKbUZKARWaq0/i8Z5hRBCdJ9OtxAAtNYh4CylVCLwplJqnNZ6e8NjlFLzgfkATqdzSm5ubjReWgjRARs3bizWWqf1dBwiNnV6DOG0Eyr1M8CttX60pWOmTp2qN2zYENXXFUK0TSm1UWs9tafjELEpGrOM0mpbBiil7MClwO7OnlcIIUT3ikaXUQbwZ6WUkUiC+ZvWekUUziuEEKIbRWOW0VZgUhRiEUII0YOiMqgshOhfNm7cmG4ymZ4HxiEVD3qLMLA9GAzeNmXKlMLmDpCEIIToMJPJ9PzAgQPHpKWllRkMBll31AuEw2FVVFSUd+LEieeBOc0dI5ldCHEmxqWlpVVKMug9DAaDTktLqyDSqmv+mG6MRwjRdxgkGfQ+tf9mLV73JSEIIXql/Px888UXXzw8JydnXFZW1ribb745y+v1KoAnn3wy5cYbb8zu6RibcjgczU7AMRqNU3Jzc/NGjBgxdvTo0XlLliwZEAqFWj3Xnj17LH/84x+ToxmfJAQhRLcoKCgwn3322aMPHTrU6bHLcDjMt771rRFz5swpLygo2H7gwIHtNTU1hrvvvjszGrE2JxAIdNWpsVqt4d27d+/ct2/fjg8//HDvypUrExYvXjyotZ/56quvrMuWLZOEIITofe67776MjRs3un7yk5+0eqFrj3feeSfOarWG77777hIAk8nEH//4x8PLli1LraqqMgAcPXrUPG3atNE5OTnj7rnnngyAyspKw4UXXjhi9OjReSNHjhz73HPPJQF8/PHHjrPPPnv02LFjx1xwwQUjCwoKzADTpk0bfcstt2SNGzduzI9+9KOMQYMGja+7c6+srDQMHDhwgs/nUzt27LB+/etfHzl27NgxU6ZMGb1p0yYbwO7duy1nnXVW7qhRo/Luuuuudr3vzMzM4PPPP3/wxRdfTA+Hw+zZs8cyZcqU0Xl5eWPy8vLGrFy50ln7+8zcsGGDKzc3N+/BBx9Mb+m4jpBZRkKILmW1Wif7/X5V9/3LL7+c9vLLL6dZLBbt8/m+PJNzbtu2zT5x4kR3w8eSk5PDGRkZ/p07d1oBtm7d6ty2bdsOl8sVnjRpUt5VV11VsX//fsvAgQMDq1ev3gdQUlJi9Pl86q677sp+99139w0aNCj43HPPJS1evDjz1VdfPQjg9/vV9u3bdwFs3rzZ8d5778XNnj27atmyZQnTp0+vsFqt+rbbbst59tlnC8aPH+/78MMPnQsWLMhev3793oULF2bfdtttRXfccUfJL3/5y3bXkMrLy/OHQiGOHj1qGjRoUPDjjz/e63A49LZt26xz584dtn379l0PP/zw0ccee2zARx99tA+gqqrK0NxxHfm9xnRCKCipwWIy4DCbcFiNmI3SoBGit9m7d++2O++8c/DKlSsTvV6vwWazhS+77LLy3/3ud4e78nUvuOCCyoEDB4YArrzyyrLVq1e7vvWtb1Xcd999WQsWLMi86qqrKmbOnFn9xRdf2L766iv7RRddNAoi3VFpaWn1/UNz584trfv6uuuuK1u6dGnS7Nmzq/72t78lL1y4sKiiosKwadMm13XXXTe87ri6BPjll1+6/v73v+cD3H777SUPPfTQ4I6+D7/fr2699dacnTt32g0GAwUFBdbOHNeamE4I/9hxghrfqYEVi8mA02LEaTXhsppw2SKf42xm4mwm4mwmHJaYfktC9Ds5OTmBuLi4kM/nM9S2CgxxcXGh7Ozs4Jmec9y4cZ7ly5cnNXystLTUcPz4cUteXp7vs88+c0T27jpFKcWECRN8X3755c7XX3894ac//WnmqlWrKr/zne+UjxgxwrN58+Zma7DFxcWF676eO3du+UMPPZR58uRJ4/bt2x2zZ8+urKysNMTFxQV37969s7mfP5PZWDt37rQYjUYyMzODixcvHpSenh54/fXXD4TDYex2+5Tmfubhhx8e0J7jWtOrbrn9wTBl7gBHyjzsPlHFhoNlrN5TxDtbjvHXzw7xzL/289SHX/HSuoMs33SUD3efZGNBGfsKqymu9hEIhdt+ESFE1BUVFZlvuOGGotWrV++64YYbigoLC82dOd+cOXOqvF6v4amnnkoBCAaDLFy4MOu6664rrruAf/LJJ/EnT540VldXq/feey9x+vTp1QcPHjTHxcWFFy5cWLpo0aITmzdvdkyYMMFbWlpqWrVqlRPA5/OpDRs22Jp73YSEhPCECRNqbr/99uyLL764wmQykZycHB48eLD/hRdeSIJIC2PdunV2gMmTJ1c/99xzyQDPPfdcSnve27Fjx0w/+MEPcm6++eZCg8FARUWFMSMjI2A0Gnn66adT6sYwEhISQtXV1ca6n2vpuI7oc7fTgZCmpNpPSbX/tOeUApfVRJLDQpLTTJLDQrLTQpLTQrytU3+fQohWfPDBB/l1X5977rmHOns+g8HA8uXL982fPz/n17/+dUY4HOaiiy6qePLJJ4/WHTNhwoSaOXPmDD9x4oTl2muvLfnGN77hfv311+N//OMfDzYYDJhMJv30008X2Gw2/corr+Tfdddd2VVVVcZQKKQWLFhwcurUqd7mXvs73/lO2S233DJsxYoVe+oeW7p06f4f/OAHOb/61a8ygsGguvrqq0vPPfdcz9NPP33o+uuvH/bb3/524MyZM8tbej8+n8+Qm5ubFwwGldFo1N/97ndLlixZchLgP/7jPwqvueaa4a+88krKRRddVGG328MA06ZN8xiNRj169Oi8efPmFbd0XEdEfT+E9mjvfgjPrslv1GXUlSwmA8lOCylOCykuK6kuC6kuK05rn8uZoh+L1n4IW7ZsOThx4sTiaMQkuteWLVtSJ06cOKS55+RqV8sfDHOiwsuJisY3BQ6LkbQ466kPl5Vkp4Wm/ZNCCNHbSUJog9sfoqDETUHJqRluFpOBNJeVtHgrA+NtDIi3keQwS5IQQvRqkhDOgD8Y5mi5h6PlnvrHrGYDA+JsDEyIfAxKsGO3GFs5ixBCxBZJCFHiC4Q5VOrmUOmplkSiw0xGgp3MRDsZiTZSpKtJCBHDJCF0oXJ3gHJ3gF3HKwGwmY0MSrSRmWhnUKKdAfE2jAZJEEKI2NDphKCUygJeAgYAGnhWa/1EZ8/bF3kDIfYX1bC/qAYAs1GRkWBncJKdwckOBkqCEEL0oGi0EILAPVrrL5VSccBGpdRKrXWzq/bEKYGQPtXNlF+C2agYlGgnK9lBdrKD9DirdDEJ0QKl1JQ5c+aUvvXWWwcgUo00PT194llnnVVTV99HdEynE4LW+jhwvPbrKqXULiATkITQQYGQbjSjyWY2MjjJTnayg5wUB4kOSw9HKETssNvt4T179tirq6uVy+XSb775ZvyAAQO6rkZ1PxDV0hVKqSHAJOCzaJ63v/IGQuwrrObD3YW8+OlB/veTA6zaeZKvTlbhDXTPgj0hYtkll1xS8eqrryYCLF26NPmaa66pL0RXWVlpuO6664aMHz9+zJgxY/L+8pe/JAItlpNesWJF3LRp00bPnDlz2NChQ8fOmTNnaDjcv8rdRG1QWSnlAl4H/kNrXdnM8/OB+QDZ2TG3kVGvUOkJsO1oBduOVmBQioEJVnJSnAxNdUr3kugxt9xC1vbtOKJ5znHjcL/wAm1WQ/3+979fumTJkozvfve75bt27XLceuutJWvXrnUB/OQnP8mYMWNG5auvvnqwuLjYOHXq1DFz5sypbKmcNMCuXbvsmzdv3j9kyJDAlClTcleuXOn65je/WR3N9xbLopIQlFJmIsngZa31G80do7V+FngWIqUrovG6/VlYa46VezlW7mVdfgkOi5GcFAdDUp0MSXFiM8saCNH3nXPOOZ4jR45Yn3vuueRLLrmkouFzq1evjv/HP/6R+OSTTw6ESNG6ffv2WXJycgItlYkeP358zfDhwwMAY8eOdefn5/erftpozDJSwP8Cu7TWv+l8SOJMuP0hdh2vYtfxqvrWw5AUJ0PTnKTHNVu4UYioaM+dfFeaOXNm+ZIlS7I++OCDPYWFhfXXNK01r7322r6JEyf6Gh6/aNGiFstJW63W+ptVo9FIMBjsV83uaIwhnA98H7hIKbW59uOKKJxXnKG61sPa/BJeXn+I5z/ez6qdJ8kvqpYS4KLPWbBgQfHixYuPTZs2zdPw8RkzZlQ+9thjA+rGAT799FM7RKdMdF8VjVlGnwD9Kov2NlXeYP3Yg9moyEp2MCzVxdA0Jy6p5ip6ueHDhwfuv//+wqaPP/LII8fmz5+fnZubmxcOh1VWVpbvo48+2heNMtF9lZS/7seUgvQ4G0NTnQxPc5IeL11LfZ2UvxZS/lo0S2s4WenlZKWX9ftLiLOZGJ7mYliak8FJDlk1LUQ/IwlB1KvyBtl8uJzNh8uxmAwMTXUyLE1mLQnRX0hCEM3yB8PsOVHFnhNVGA2KzEQ7w9MjrQfZblSIvkkSgmhTKHyq5tJHuyE93lrftSRTWoXoOyQhiA4rrPRRWOljXX4J8XYzw9OcDE9zkZloxyDjDkL0WpIQRKdUegJsOlTOpkPl2MzG+hlLOSlOLKaolsoSQnQx+R/bgypLCnnqnu9RWVrU06FEhTcQYtfxSlZsPc4z/8pn+aajbDtSQY0v2NOhiT7q0KFDplmzZg3LysoaN3bs2DHTp08fsXXrVmvbP9nY+++/7xoxYsTY3NzcvAMHDphnzpw5rLnjpk2bNnrNmjVRrdsUS6SF0IM+ePlpDmzfwAd/+T3X3vVAT4cTVcGw5kBxDQeKa1C7YUC8jWGpToanu0h1dfj/q4hxD63YOSia5/vprLxjbR0TDoeZM2fOiHnz5pWsWLFiP8C6devsx44dM0+YMMHX1s839NJLLyUvWrTo+MKFC0sB3n///f1nFnnvJi2EHnDvrAksumw0a1csRWvN2hVLWXTZaO6dNaGnQ+sSWsOJikgpjf9bV8ALnxzgoz2FHCpxEwpLnUNxZlasWBFnMpn0vffeW9/EPvfccz2XXXZZ9e233z545MiRY0eNGpX33HPPJdUd31x569/85jep7777bvLDDz+cOWfOnKF79uyxjBw5cixAdXW1mjVr1rBhw4aNvfTSS4d7vd76QbI33ngj/qyzzsrNy8sbc/nllw+rqKgwAGRmZo7/z//8z0F5eXljRo0albdp0yYbQEVFheHaa68dMmrUqLxRo0bl/elPf0ps7Tw9QRJCD7j/z6uYPGMWZmtkho7ZamPyRbO5/6V/9nBk3aPCE2DzoXJe//IIz6zJ592tx9l5rBKPX1ali/bbunWrfeLEie6mj7/00kuJ27Zts+/atWvHP//5z70/+9nPBhcUFJghUt7697///eF9+/btOHTokHXlypWuRYsWFV9yySXlP//5z4+8/fbbBxqe69FHH0232+3h/fv37/j5z39+bOfOnU6A48ePm37xi19krFmzZu/OnTt3TZ482f3QQw8NqPu51NTU4M6dO3fdcsstRY888sgAgB/96EcZ8fHxob179+7cu3fvziuvvLKqrfN0N+ky6gHxKenYHC6Cfh8mi5Wg34fN4SI+Oa2nQ+t2vkCYvSer2HsyUqU1I8HG0NrFcGlx0rUkOu7jjz+O+853vlNqMpnIysoKnnPOOdWffPKJIyEhIdzR8taffPKJ66677iqESKntUaNGuQFWr17tzM/Pt02bNi0XIBAIqClTptTvmzBv3rwygGnTprnffvvtJIA1a9bEv/LKK/VdUWlpaaGlS5cmtHae7iYJoRtUlhTy0i8WcfXC+3jz6Ye58b7HqSov4bxZc/naFd9l/XvL+szAcmeEteZouYej5R4++aqYOJupfqV0VrIDs1EatOKU8ePHe5YvX57UkZ+JVnlrrTUXXHBB5TvvvHOguedtNpsGMJlMurXXaOs83U3+h0VRS7OG6gaP//LI4vpB5JuXPMU1dy4hc3gu19y5hJuXPNVDUceuKm+QLYcreGvzMf64Op83vjzCxoIySmv8PR2aiAGzZ8+u8vv96tFHH02te+yzzz6zJyYmBl977bXkYDDIsWPHTJ9//rnr61//es2ZvMYFF1xQ/fLLLycDfPHFF7a9e/c6AC688MKaDRs2uLZv326FyHadbc1umj59euXjjz+eXvd9UVGR8UzO05UkIURRw1lDcPrg8cmCfR0eRO5rU1PPVDCsKShxs2ZvEX9ee7B+j4d9hbK/dH9lMBh4++238z/88MP4rKyscSNGjBj7wx/+MPOmm24qHTt2rGfMmDFjL7zwwlEPPvjgkezs7DOa+7x48eLCmpoa47Bhw8bed999mXl5eTUAgwYNCj7zzDMHr7/++mGjRo3Kmzp1au62bdtaXbb/y1/+8nh5eblx5MiRY0ePHp333nvvxZ3JebqSlL+OgntnTSDoP32Wm8Fk5qyvf5Otn65s9LzJYmXCBZcxZ/4P2xw3eO3JB1j37iuce+X1fW5qarTU7RCXlewgJ8VJRrxNVky3QMpfCyl/3cXu//Mq3n72V2xbu4qAz4syGNHhEJNnXInFYiMU8KMMBnQ4jDIYCAX8bQ4i33vleIKBU10ja1csZe2KpZgsVv5nxdbueFu9RsP9pT/bX4rFZCAz0U5WsoOsZDtpLiuRnV6FEK2RhBAFdbOGAj4vADocadVsWLkcAGUwMHzCNKrLSnAlpTAga1ibXUBnTb+CDauWYzAYCYdDmK02xp9/KTOuu5Wn7vkeN973eL+cldQe/mC4flEcgN1iJDPRzuAkO4OTHKS6LJIghGiGJIQoqSov4exLr6aqvIS9X35KOBTCZLFiszu5/Vcvkjkst13nadr9FK5NLgGfF5vDxbp3l/XZ1c1dxeMPsa+wmn2Fkdl8douRQYl2Mms/0uOs0sUkBFEaVFZKvaCUKlRKbY/G+WJV3QDv0fxdpw303rzkKeb+1yMkpw9Ch8P16wuqK0pZt+KVdp/77ieWNVq0pgxGcs/+BspgqO826g+rm7uSxx8iv7CaNXuLWPr5If7wr3xe3XCYT/cVs7+oWgap2yccDocli/Yytf9mLe4hHa0Wwp+Ap4CXonS+mFK3jiB5YGb99NHCQ/nN3qVXlZeAUo3u8tvT/183Q2ndu8saLVoLBfwkD8hkycv/ajROUdeFNGf+D7vyrfcL/mCYI2UejpR56h9LcpgZmGBjQLyNjAQ7qS4LJlkH0dD2oqKivLS0tAqDwSD1R3qBcDisioqKEoAWb9yjkhC01muUUkOica5Y9OAN09HhMPu3fQHAyYJ9QPMX+puXPEVlSWG7L95Nu4jWrlga+UIZuPuJv9UvWmtpdTNaNxpTqEteMsbQOWXuAGXuALuOVwFgNCiSnRYGxNtIj7OSFmcl1WXttyW+g8HgbSdOnHj+xIkT45Dp671FGNgeDAZva+mAbhtDUErNB+YDZGdnd9fLnpG6i+qh3VsazfRpquH00YY6UpqipRlKUy+ZU79orU5zq5ubVkztyxVUe1IorCmq8lFUdSp5KwUJdjNpcVZSnFZSXRZSXFYS7eY+PyYxZcqUQmBOT8choqvbEoLW+lngWYisQ+iu1z0TdRfVKRdfRTgUrL9YN9TW9NH2lqZobYbShpXLT2t91Pns/deanZba9HuZptp1tIZyd4Byd4CvOFV+xmRQJDrMJDktJDkiH4kOM4kOMw6LzOMQsatf/HW2txulaffNhlXL679WyoDWYRJS0rG7EtqcPtrw4t3wLr85zc1QamuMoOm0VJPFSlxSKlVlxQT9Phlj6EHBsKa42k9x9emtS4vJQILdTLzdTLzNVPvZTJzNhNNqwmkxypRY0WP6RUJobzdKS903aYOHcON9T9Tf6Ue77lDd+V57YkmjGUoNWx8tdWPVTUsN+n1YrHZCAX+/r6Aay/zB8GldTw0ZlMJpNeKwmHBYjNgtRpwWEzazAZvZiM1swGoyYjUZsNR+mI0GKfwnoiIqCUEptRS4EEhVSh0Blmit/zca5+6M5gZsW+tGaan7pujIQR5bcFWXd7+01s3UUjeWMhgZPeV8XIkp7N34iVRQ7eXCWlPlDVLl7Vjpnf+4ZKS0LESnxWwto9/9Dv76bhVmexCbM4zdGcbuCuGMD+GID+GMD9d/bWjh5qi12T4t3Tm/+OAd2ByuFrtvuvuOu6U6SUD9tFSpcyTamxCiVctI9E0x22WUnw97N9twVxnwuo0tHmcwaJyJIeKTgsQlhYhPCZKQEiQhNUhCipNQMI+A7xNMFt2ubpT2dN/U6Y4pnk27scxWG464REZOOpfp375JWgJCiKiJ2YTw299C3rcPUeMLEQ6Bz2PAXWXAXWWkprL2o8JIdbmRqjITVWWRz8f2W6kqN6LrF1H+CvgVBkMAi+sE29YeJhx28tWmvzLzxivIznWQPDCA0dj4At+eWUJ13Tgrnn+U0pNHuyQxNDeFdezXZtS3CNoasBZCiPaK2S4jOLPy15Ulhfz54f/i6oVPEgoMoKzIRHmhmdKTJspOmik5YabwsCYcOlVy3GjSpAwMEAhso+zkPxkxMZ4rb72CAdk+bI7Tfz8tdeMog4HH3t/VoXjb48UH7yA+Oa1RcpINdURD0mUkoqHPJYTW9g9ofCEfAAwHRgKjaj9G134+tWFR0oAAGUN8ZAz1kTHUz6BhPmyOI6x4/ld8+dGKZmOIhbn/smK5f5GEIKIhZruMOqo9M4oa98efRBmK0eG1TPjGTEzGdbWPBzCZc8nOncvQsTdQejKZ4wes7N7gJByK/IczW7Ox2h8HrgQ2AZ8DOzBbjTEz9787urOEEH1Ln0kIzQ2+Nndx/mrLZ6dNK9265v365yMzd3YyMGcHV97iB04AEAxA4WELR/OtHMu3sWFVNUbz9wkF/q32J30EfFs4ll/O3i+HkZ3rJS0zQHM3bV1x997SOoW6xXUPzvsGj72/S1oOQogW9ZmE0J76QR+8/DRVpUWkZg4hdVB2/bTSlhagNWQyw6BhfgYN88OlVVz1bw60PsozP34YmEp8yuXs/dJF4eFx/PV/HADY40Lk5HoZmudhSJ6H7FwvVrvuknpDTdcpNO3O0uEwiy4bjTIYQOv615YEIYSo06fGEFoafG1tLn9zOjMGEA7BiUMWDu22UbDbTsFOGycK6sYkQsBmYA3wMfAJUNSp12vtvdWV22iVMqDQspahl5MxBBENMZsQjh8/zvTLr+KGH3f+zrW5BWqOuERy8ibh99R0+QI0d5WBgt029mwI8+VH5VSXjwTsANgchxh3npkx0zTDJ3iIT+74rKqW1in4vR62fvx+fWJIHZRDRcnJ0wr1NRQLA+Ki4yQhiGiI2QIoDz30EPu2RrpVGqrbWawji7Famst/0/2/PW2Hs66o/+OICzPmbDffWuBlwvm/BBIxmqYDP8LmKmHb2nT+7xeDeOD64Txyaw6vP5XO1k9c1FS2/c/T0nub91+PoMMhzp89j0VPv8n5s+cRDgdbTQaTL5rN/S/9M3pvXAjRq8TcGILdbsfrPXXRajpb6Ez731taaNbeMtXRUlVewvmzr23wej/kxvuf4ug+K/u2ONi3xcEXH8Tz6duJKKXJHO5j5GQ3oye7GTrOg9lyeouupffQtOLqiw/ewfDx0xqV5YDI+gkdDrNv8/oufe9CiNgWc11Gx48fZ/HixSxfvhy3213fjbP14380u1mN0WwhJ3dinxoUDQXh0B4bX21ysHeTg4JddkJBhckSZtg4D7lT3YyaUkPGEH+zs5ia03Dw+IP/e4p17y0DFFqHmfiNmZws2MeJgn2cN2uujCX0QtJlJKIh5loIGRkZxMfH4/V6MTfoxrn/pX82O63UYDSycdVbfWqHMKMJho71MnSsl8u+V4rPo8jf5mDvRgd7vnTw9rNpQBrxKUFyp9aQO7WGUZPdOOJaHkBu2LJq2KL4zb9fzZYG025lYx0h+q+YayEAfPvb3yYjI4PkKZfz0fKl9bOFXntiCeveW4bRbGm1Amhfv5CVFZrYs9HBno1O9n7pwFNtRBk0Q8Z4yZ1aw5hzasgc7kOplmch1f2eWqoIO+O6W3nz6Yf7VMurL5MWgoiGmGshALzxxhtAZNrppfMW8NIvFlFZWtToznbV0j+yfe0qlMHQ73YIS0oP8rXLK/na5ZWEQnBot43dXzjZ9YWTv/85lb//OZX45CBjptVwzR2b2P3FEnZ+/k6zC/aaDkoHfF72bV6PwWiUvZmF6GdiMiE01LAEQ01lOdfcuYT45DSc8YmEggGAfr1DmNF4qnvp8ptKqCozsvsLJzs/d7LlYxefvT8GZXgZHf4Ig/F9Ar63Tvs9NUy0jy38FpWlRWxYuRyQLiQh+pOY7DKC02cbtUUZDIw792KpAtpAKAgHdtj52+Prqak8D091JgBW+0HOnxPP2K9Vk5PrxWBsfYHb5Itm98jmQKL9pMtIREPMthD279/Pt268nc9XvdPqcT25m1msM5pgxEQPP/nTRKCGoqMH2PmZk53r01n9mp0PlyXjTAiSN62Gb//7JvZs+Bk7P19Rv1ahrsxFf2x5CdEfRWtP5ZnAE4AReF5r/Uhnz5mRkYHN4ao9/6kSDGarjYSUAZQcP1Q/uGxzuEBrnrrnezII2oq0zADTv13O9G+X46kxsGeDg+3rXGxf7+KLlXVdSyuBFcBbTLhgIq6EZNmRTYh+otMrlZVSRuD3wOVAHjBXKZXX2fMCVJaVcP7seYy/4LLa14oMIIfDIc6bNZe7n/gb582aS1VZcaNplaJtdmeYs6ZX870fneC/l+Wz8NeHSR74DjbHVOAPwDH2fvk48SlPMPPGZ2nYs3gmq8WFELGv02MISqlzgQe01t+s/f7HAFrrX7b0Mx0tbtfajmFtTasUHaM1nCywRFoOa50c2hOpuZSS4WfcuTWMPa+aTR/+mPV//2t9QTypmNrzZAxBREM0EsK1wEyt9W21338fOEdrfUeT4+YD8wGys7OnFBQUtHnu9lQ7bWkevYwpREdliZEd611sX+dk1+dmIrvJFRPpVloOrEQZvKClYmpPkoQgoqHbittprZ/VWk/VWk9NS4vehbo9+yCIzjjOxg9n891FW7jzNyuxOm5CGT4AriKSEIrR4TfR+ibWrljFostGc++sCT0bshDijERjUPkokNXg+8G1j3Wb7i5Q1580HZvxuZcCf8ZodhIKTMPm/D7emouBOUAYg/FzzpsVR+ERM+mDAz0Zer9RWVLIhRf+gGXLljFw4MCeDkf0YtHoMjIBe4GLiSSCL4B5WusdLf1MV22Q05D0a3dOezYVUgYDSekZlJ44htE0lVDwcuBbwCQA0rN8jDuvhnHnVpOd68UQs8XWe7fXnnyA9e8t4/bbb+fpp59u9VjpMhKticrCNKXUFcBviUw7fUFr/XBrx3dHQnjtyQdY9+4r0q99hpqOzSiDAVDocKhRvaM//PAmPNWV6HDDwnrZwBxQV2MwzCAcUrgSg4z9Wg1jz61m1CQ3Flv3L4js7epucq5eeB9vPv0wBbs216/Wb8hms+HxeJo9hyQE0ZqorEPQWr8HvBeNc3VW0ztbKb1wZpobm4HGZULWvbsMT1VF/T7Op5LHUXT4KaZecoSrF2Sz6wsnO9Y52bLGxWfvJ2CyhBk1yc3Yc2sYe0418SlnlvT7iva2Zuu67/7yyGIKD+WTMiiH4qMHMRiMhMMhHA4HV199NY8++mg3Ri/6kpgtXQFn1kKQWUfR03C674sP/jsANy/5Pb/596ubtAhaV5eMgwHYv83BjvVOdqxzUXrSDEDWKC9551Qz9ms1ZI7wtXuPh76irdZsR/YEX7BgQavdRtJCEK3pcwkBaFQmOxTwS7dRlEVjj2qt4fhBCzvXu9j5mZOCXTa0VsSnBMmbVk3eOTWMnOTGau+7XUstXeibbvpU9/ve+unKZo9XBgMzv/lN0tPTqaysrK8W3BxJCKI1MVvLqDNk1lHXamkf52vveoDXnljSrj2qlYJBQ/0MGlrKJXNL66u07ljvZNO/4lj/90SM5jAjJngYM62GMdNqSMvsG7OW6rqIbv3vP/CXXy7G56lpVMK96aZPdb/vUMBfv91pnbp6Uzk5OfzhD3/owXcl+oI+2UIQXa+l1eOtrSqv01afeTAAB7bb2fm5k12fuSg8YgEgdZCf3LNrGHN2DcMneGJ6YLq191jXRZSePZyTBfsAGo3TNGWyWBlz9jeIT07j5OH9VJeVUFVeTFxiKq6kFAZkDSPd5Gm1ZVBHWgiiNZIQRJdp6aLY0RlgJcfN7PrCwa7Pnezb4iDgM2Ayhxk23lO7haibATnt31+6OzT3HtszFmA0mZvd9Kmt8S9ZqSyiQRKC6DJNL4rRqDvl9ykObLOze4OT3RsdnCywAhCfEmT05BpGTXEzapKbuKSe+btpbVzgp7X7gjcdCzBZrEy44DIMRmP9xkQmi7V+/OuyGxa2OQtJEoKIhj45hiB6VktTfw0mM5NnzGp2Blh7Waya0VPdjJ7q5ipO7S+990snO9a7+GJlAgAZQ32MmuRm5CQ3w8a7sTm658bn/j+varJ+w4gOh5h04RXNjgXUtQa+/LDxvh9Bvw9lMJxWyVcmR4iuJAlBRF1LF8XJM67EYrFFte5UUnqQvGn72LBqEYufeZzKksHs/dLBV5sdfPpOAv96IwmDQZM12svwCR5GTnQzZKyny2Yv1V306zYZ0uFIS2XDyuVsWLkcpQycN2tu/ViAKymFpLQM9m5ah7uqvFGi3PrxP9j26cr6c8uaGtHVJCGIqGvtogiRmTF3P/G3+kHnpitwO1pupO4OetVfI3fQ2aN9XDK3jIBfcXCHja82O9i3xcHq15L4cFkyBqMma5SXYeM9DB/vYeg4D3Zn+9dVtKWqvISzL72aqvKSdk3BhVNTpRsmyvtru5g606ISoiMkIYgu0Z6L4qXzFvDSLxax4oXHGq3AbW/XSFur0s0WzchJHkZO8gAl+DyKAzvs7NviYP82O2veSOKjvyWjlGbgUD9D8yLJYehYD0npwTMepK6bVdXeKbh1v6+mU6Wlkq/objKoLLpUw0WCQb+P+OQ0Fj39JvHJadwzc0yrK56bdo00nbXU2VXpfq+iYLeN/dvsHNhhp2CXHZ8nUoEvPjlIzhgPOWO8DBnjJXOEt8PdTO2Zghutc8igsogGaSGILtXwzvflR+7hRME+Hpz3jTYTwYQLLmPGdbfy1D3fq+9KSh6Y2WhwtbN30BabZuRZHkaeFSkEFw7B8YNWDuywUbDLTsEuG9s+jQNAGTQDsv1kj/aSNcrL4JE+Mob6sFhbThINL9zX3LmkXTF1xTmEaC9pIYgu15FaPHUrb8+98noA1r37Chqgmb/Thgu2OnMX3prqciMFu20c3mvj0B4bh/fYqKk01j4bJD3bR/aoIIOG+Rg03MegYT5cCdEbj2gvaSGIaJCEILpcc107CSkDKDl+CFBoHSYhJR27KwFXUgr5Wz9vs3je5Itm90jBQq2h9ISJ1558iz0b3SSmXUkoOI6qslON7fiUIANzfAzM8TNwSOSz3XWUV5+4u0MD5h3Z00MSgogG6TISXa65rp1wOHTaIGrdnX1rxdzqWhBnOrja2Y2Tfji7cWunvOh+AIzmLH7w0FqO7bdybL+FEwVW1r2XQMBXtytQNvAOv1lYSt456aQNDpCSESA1009qRqDZMhyy/kB0N0kIols0N4umrk+8ad94a8XcJlxwGa6E5EbTVZu7uLf0XGcvsk3XWDQeyHYzMOcg29ct4raHHseVmMYPZ19FKDACyAVGU1may/q/pwGN441LDpI8IEDygACb1zyDDu8HqoCJrF2xUtYfiG4hCUF0i44OjtYlkIYLuAZkDWuUSF578oEWL+5NL/zR2jiprYHspq/705eeq00gTzVKIJfe8BP8ngyKj5kpPmah9KSZshMmDu21odQ9aBrvN6pUAEdcmN/eGSY+OYQrMYgrMYQzIYQrIcRaB5x/frvfhhDNkjEE0evce+V4ggH/aY+bLFbQutnnjGYLEy+4LCobJzU3FXTXF2tarNM07dKrO7Q/RzgErzz6OBv+uRGDaSjhYAaZwy8iY+h0KktNVJaaqC43UlNpRIcj4wbTp2tWr5YxBNE5nWohKKWuAx4AxgDTtNZtX+WF6KSzpl/BhlXL67eObHhxX/H8o80+N+O6W3n2vh/U39kHfF72bV5/Rq/fXGuntTURr//uwQ7tz2Ewgs+bz/mzRzb4maeYd29uo+PCYfBUG6ipMHLjuUPO6L0I0VCnWghKqTFAGHgGWNzehNDeFsLa/GI8/hDBsCYU1gRCYQKhus9h/MEwvmDks+j7Wpu+2nSsoaHzZs0FIt1EA3JG8L0fPVa/JuK8WXO59q4HOj3YDD27U5/MMhLREJUuI6XUarogIbRXKKzxBUN4A2Hc/iAefwhPIESNL0SNL0iNP0iNL0S1L4DbH2puSrvoBZrehSuDkdFTzseVmEJ1eQkOV/xpz+3Z+Em79n9uuP7hTC7ilSWF/M/ts8k7ZwbTv31Tl6yJaI0kBBEN3TaorJSaD8wHyM7Ojuq5jQaFw2LCYYFkp6XVY0NhTZU3QJU3SIUnQKUnQEXtR7kngMcvYxaxqumAbijgJ3lAZv0F/LUnlpz23JKX/9UoiZgsVuKSUqkqK27U2qhLGk0Hm5srl9HS7CVPVQUWq43M4bmyqlj0Sm0mBKXUKmBgM0/dp7V+q70vpLV+FngWIi2EdkcYZUaDItFhIdFhIauZ572BEGVuP2U1Acrcfkpq/JRW+yj3BKRlEQNa2y+7PQXiQgE/FqudUMDfaNvK5sYj4NSsoRXPP0rpyaOnlc9oOsAtJapFb9Ynuoy6QzAUptTtp7jKT3G1j+JqH0VVPtzSooi6aPTnN9R0VtC2T1dSVV7SYldSfHIa7qqKZmcrnUYp0LrZhNKdq6ily0hEg6xDaCeT0UB6nI30OFujx6u8AQqrfBRW+iis8nKy0itTZTspmit0K0sKqaks55o7lxCfnMY1dy7hmjuXtDgeUXL8MEVHDjDh6zPJ3/o5NRWlrb9A7Q1VuHbPh4DPKyWqRa/V2WmnVwO/I7Ls8l2l1Gat9TejElkvEWczE2czMzzNVf9YpTfAyQovJyq9HK/wUljpJRCS/qa2RGvxWEMtJZfmFpjt/mJN/fNbP36//mulDGh9qjVRN6MpdVAOFSUnTxvgriorPqNYhehpnUoIWus3gTejFEufEW8zE28zM3JApHRyOKwpqvZxvMLLsXIPx8o9VHmDPRxl7GmtLERHtSe5NBxzWPPGn9i7aR0VxSdOO1fDZACnymfs+uJfLQ5wC9EbSZdRNzAYFAPibQyIt3FWViIAFZ4Ax8o9HC3zcKTMTZk70LNBxoBo7hDWnuTScEro3P96hNeeWMLad19p1CIwW2044hIZOencRtNJr7lzCS8+eAdjzp7e7gVnQsQ6SQg9JMFuJsFuZkxGPADVviBHyzwcLnVzuMxNeT9NEK3NIuqIM0kuVeUlnD97HlXlpWz9+H2UMhD0+xj7tRn1d/4Np5PK5jWir5GEECNcVhOjB8YxemCkm6nSG4gkh9JIkqj29Y8upmheZDuaXOpe+8UH7+D82fPkzl/0OzFd3E6cUlzt41Cpm8Olbo6UeaRch2hEpp2KaJAWQi+R6rKS6rIyOTuJUFhzrNzDwZIaCkrcFFW1b3tKIYRojSSEXshoUGQlO8hKdvD1kZHxh4PFkeRQUFqDLyCtByFEx0lC6ANcVhPjMhMYl5lAOKw5VuHhQHENB4trKK5ux2pbIYRAEkKfYzAoBic5GJzk4Osj06j0BjhQVMOB4hoOl7oJhmWBnBCieZIQ+rh4m5mJWYlMzEokEApzqNRdnyD6y8wlIUT7SELoR8xGA8PTXAxPc6G15mSlj/3F1ewvqpGBaSGEJIT+SinFwAQbAxNsnDc8lQpPgP1FkeRwtNxDSLqWhOh3JCEIILJyelJ2EpOyk/AGQhSUuMkvquZgicxaEqK/kIQgTmMzG+tXTYfCmiNlbvYX1ZBfVC1F+YTowyQhiFYZDYqcFCc5KU5m5KZTWOllX1E1+UU1FMu4gxB9iiQE0SHp8TbS42vHHdwB8ouryS+s5li5l7DsMSpEryYJQZyxBIeZydlJTK4dd6jrVjpU6pZaS0L0QpIQRFTYzEbyBsWTNyieYO16h/1FNewvrpYtRYXoJSQhiKgzGQ0MS3MxLM2F1umcqPRGkkNRtZTSECKGdXZP5V8DswE/kA/crLUuj0Jcoo9QSpGRYCcjwc75I1Ipd/vZX1wTWe9Q5pFxByFiSGdbCCuBH2utg0qpXwE/Bjq+Aa7oNxIdFiZnW+rHHQpK3OwvquZgiRtvQLqWhOhJnUoIWusPGny7Hri2c+GI/qTheoeGVVoPFNdQIl1LQnS7aI4h3AIsi+L5RD/StEprhSdQX8JbqrQK0T3aTAhKqVXAwGaeuk9r/VbtMfcBQeDlVs4zH5gPkJ2dfUbBiv4jwW7mrKxEzqqt0nqkzMPB2tZDhSfQ0+EJ0Se1mRC01pe09rxS6iZgFnCxbmWDZq31s8CzENlTuWNhiv7MbDQwNNXJ0FQnM4CyGn/99qFHytwEQvLnJEQ0dHaW0UzgXmC61todnZCEaF2S00KS08Kk7CSCoTDHyr0UlEYSRHG1D5m4JMSZ6ewYwlOAFViplAJYr7X+t05HJUQ7mYwGslMcZKdE9peu8QU5VOrmUKmbw6VuKcYnRAd0dpbRiGgFIkQ0OK0mxmTEMyYjHoDSGj+HS90cLnNzpMyDxy9TW4VoiaxUFn1astNCstPCxKxEtNYUV/s5XObmaJmHo+WSIIRoSBKC6DeUUqTFWUmLszI5OwmtNSU1fo6UeThWHvmQLibRn0lCEP2WUopUl5VUl5WzshIBqHAHOFru4XiFh2MVXkpkkFr0I5IQhGggwWEmwWEmb1BkDMIXDHGywsfxCg8nKr2cqPDilm4m0UdJQhCiFVaTsX4WU50Kd4CTVV5OVno5WenjZKVX9n8QfYIkBCE6qK4VMWpAHABaayo8AQqrfBRV+Sis8lJc5afaJ+MRoneRhCBEJymlSHRYSHRY6pMEgNsfpLjKT1G1j5JqHyU1fkpr/NKaEDFLEoIQXcRhMZGdYmrU3aS1ptIbpLQ2OZTV+Cl1Rz7L2IToaZIQhOhGSikS7GYS7GaGpjobPecLhih3Byhz+6lwByj3BKjwBKj0BKj2BWW2k+hykhCEiBFWk5EB8UYGxNtOey4U1lR5A1R5g1R4Ip+rfUGqfZGva0vHCNEpkhCE6AWMhlPjFFk9HYzosww9HYAQQojYIAlBCCEEIAlBCCFELUkIQgghAEkIQgghaklCEEIIAUhCEEIIUatTCUEp9ZBSaqtSarNS6gOl1KBoBSaEEKJ7dbaF8Gut9QSt9VnACuBnnQ9JCCFET+hUQtBaVzb41glItRUhhOilOl26Qin1MHAjUAHM6HREQggheoTSbZRQVEqtAgY289R9Wuu3Ghz3Y8CmtV7SwnnmA/Nrvx0N7GlHfKlAcTuO60mxHmOsxwexH2OsxwftjzFHa53W1cGI3qnNhNDuEymVDbyntR4XlRNGzrlBaz01WufrCrEeY6zHB7EfY6zHB70jRhH7OjvLaGSDb68CdncuHCGEED2ls2MIjyilRgNhoAD4t86HJIQQoid0KiFora+JViAteLaLzx8NsR5jrMcHsR9jrMcHvSNGEeOiNoYghBCid5PSFUIIIYAYSQhKqZlKqT1KqX1KqR8187xVKbWs9vnPlFJDYiy+RUqpnbVlPP6plMrpzvjaE2OD465RSmmlVLfPSGlPjEqp79T+Lncopf4aS/EppbKVUh8ppTbV/ltf0c3xvaCUKlRKbW/heaWUerI2/q1KqcndGZ/oA7TWPfoBGIF8YBhgAbYAeU2OWQj8sfbr64FlMRbfDMBR+/WC7oyvvTHWHhcHrAHWA1NjLUZgJLAJSKr9Pj3G4nsWWFD7dR5wsJt/h98AJgPbW3j+CuDvgAK+BnzWnfHJR+//iIUWwjRgn9Z6v9baD7xCZAprQ1cBf679+jXgYqWUipX4tNYfaa3dtd+uBwZ3U2ztjrHWQ8CvAG93BlerPTH+APi91roMQGtdGGPxaSC+9usE4Fg3xofWeg1Q2sohVwEv6Yj1QKJSKqN7ohN9QSwkhEzgcIPvj9Q+1uwxWusgkTIZKd0SXfvia+hWIndp3anNGGu7D7K01u92Z2ANtOf3OAoYpZT6VCm1Xik1s9uia198DwDfU0odAd4D7uye0Nqto3+rQjTS6VpG4hSl1PeAqcD0no6lIaWUAfgNcFMPh9IWE5FuowuJtLLWKKXGa63LezKoBuYCf9JaP6aUOhf4P6XUOK11uKcDEyIaYqGFcBTIavD94NrHmj1GKWUi0lwv6Zbo2hcfSqlLgPuAOVprXzfFVqetGOOAccBqpdRBIv3Lb3fzwHJ7fo9HgLe11gGt9QFgL5EEESvx3Qr8DUBrvQ6wEakhFCva9bcqREtiISF8AYxUSg1VSlmIDBq/3eSYt4H/V/v1tcCHWuvuWkDRZnxKqUnAM0SSQXf2e7crRq11hdY6VWs9RGs9hMg4xxyt9YZYibHWciKtA5RSqUS6kPbHUHyHgItr4xtDJCEUdVN87fE2cGPtbKOvARVa6+M9HZToPXq8y0hrHVRK3QH8g8hMjxe01juUUv8NbNBavw38L5Hm+T4ig2rXx1h8vwZcwKu1Y92HtNZzYizGHtXOGP8BXKaU2gmEgP/SWndLS7Cd8d0DPKeU+k8iA8w3deONCUqppUQSZmrtOMYSwFwb/x+JjGtcAewD3MDN3RWb6BtkpbIQQgggNrqMhBBCxABJCEIIIQBJCEIIIWpJQhBCCAFIQhBCCFFLEoIQQghAEoIQQohakhCEEEIA8P8BkmcqolPCIjQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    untrained_pred_dist = likelihood(model(test_x))\n",
    "    predictive_mean = untrained_pred_dist.mean\n",
    "    lower, upper = untrained_pred_dist.confidence_region()\n",
    "    \n",
    "f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "# Plot training data as black stars\n",
    "ax.plot(train_x, train_y, 'k*')\n",
    "# Plot predictive means as blue line\n",
    "ax.plot(test_x, predictive_mean, 'b')\n",
    "# Shade between the lower and upper confidence bounds\n",
    "ax.fill_between(test_x, lower, upper, alpha=0.5)\n",
    "ax.set_ylim([-3, 3])\n",
    "ax.legend(['Observed Data', 'Mean', 'Confidence'], bbox_to_anchor=(1.6,1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually, this does not look like a good fit. Now, let us train the model hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/50 - Loss: 0.944   lengthscale: 0.693   noise: 0.693\n",
      "Iter 2/50 - Loss: 0.913   lengthscale: 0.644   noise: 0.644\n",
      "Iter 3/50 - Loss: 0.879   lengthscale: 0.598   noise: 0.598\n",
      "Iter 4/50 - Loss: 0.841   lengthscale: 0.555   noise: 0.554\n",
      "Iter 5/50 - Loss: 0.798   lengthscale: 0.514   noise: 0.513\n",
      "Iter 6/50 - Loss: 0.750   lengthscale: 0.475   noise: 0.474\n",
      "Iter 7/50 - Loss: 0.698   lengthscale: 0.439   noise: 0.437\n",
      "Iter 8/50 - Loss: 0.645   lengthscale: 0.404   noise: 0.402\n",
      "Iter 9/50 - Loss: 0.594   lengthscale: 0.372   noise: 0.369\n",
      "Iter 10/50 - Loss: 0.547   lengthscale: 0.342   noise: 0.339\n",
      "Iter 11/50 - Loss: 0.505   lengthscale: 0.315   noise: 0.310\n",
      "Iter 12/50 - Loss: 0.466   lengthscale: 0.291   noise: 0.284\n",
      "Iter 13/50 - Loss: 0.430   lengthscale: 0.271   noise: 0.259\n",
      "Iter 14/50 - Loss: 0.395   lengthscale: 0.254   noise: 0.236\n",
      "Iter 15/50 - Loss: 0.361   lengthscale: 0.240   noise: 0.215\n",
      "Iter 16/50 - Loss: 0.327   lengthscale: 0.229   noise: 0.196\n",
      "Iter 17/50 - Loss: 0.294   lengthscale: 0.220   noise: 0.178\n",
      "Iter 18/50 - Loss: 0.261   lengthscale: 0.214   noise: 0.162\n",
      "Iter 19/50 - Loss: 0.228   lengthscale: 0.209   noise: 0.147\n",
      "Iter 20/50 - Loss: 0.196   lengthscale: 0.206   noise: 0.134\n",
      "Iter 21/50 - Loss: 0.164   lengthscale: 0.205   noise: 0.122\n",
      "Iter 22/50 - Loss: 0.133   lengthscale: 0.204   noise: 0.110\n",
      "Iter 23/50 - Loss: 0.104   lengthscale: 0.205   noise: 0.100\n",
      "Iter 24/50 - Loss: 0.076   lengthscale: 0.207   noise: 0.091\n",
      "Iter 25/50 - Loss: 0.049   lengthscale: 0.210   noise: 0.083\n",
      "Iter 26/50 - Loss: 0.025   lengthscale: 0.214   noise: 0.076\n",
      "Iter 27/50 - Loss: 0.003   lengthscale: 0.218   noise: 0.069\n",
      "Iter 28/50 - Loss: -0.017   lengthscale: 0.223   noise: 0.063\n",
      "Iter 29/50 - Loss: -0.033   lengthscale: 0.228   noise: 0.058\n",
      "Iter 30/50 - Loss: -0.047   lengthscale: 0.233   noise: 0.053\n",
      "Iter 31/50 - Loss: -0.057   lengthscale: 0.239   noise: 0.049\n",
      "Iter 32/50 - Loss: -0.065   lengthscale: 0.243   noise: 0.045\n",
      "Iter 33/50 - Loss: -0.069   lengthscale: 0.248   noise: 0.042\n",
      "Iter 34/50 - Loss: -0.071   lengthscale: 0.252   noise: 0.039\n",
      "Iter 35/50 - Loss: -0.071   lengthscale: 0.255   noise: 0.037\n",
      "Iter 36/50 - Loss: -0.069   lengthscale: 0.257   noise: 0.035\n",
      "Iter 37/50 - Loss: -0.066   lengthscale: 0.257   noise: 0.033\n",
      "Iter 38/50 - Loss: -0.062   lengthscale: 0.257   noise: 0.032\n",
      "Iter 39/50 - Loss: -0.058   lengthscale: 0.255   noise: 0.031\n",
      "Iter 40/50 - Loss: -0.055   lengthscale: 0.252   noise: 0.030\n",
      "Iter 41/50 - Loss: -0.054   lengthscale: 0.248   noise: 0.029\n",
      "Iter 42/50 - Loss: -0.053   lengthscale: 0.243   noise: 0.029\n",
      "Iter 43/50 - Loss: -0.053   lengthscale: 0.237   noise: 0.029\n",
      "Iter 44/50 - Loss: -0.055   lengthscale: 0.230   noise: 0.029\n",
      "Iter 45/50 - Loss: -0.057   lengthscale: 0.223   noise: 0.029\n",
      "Iter 46/50 - Loss: -0.060   lengthscale: 0.216   noise: 0.029\n",
      "Iter 47/50 - Loss: -0.063   lengthscale: 0.208   noise: 0.030\n",
      "Iter 48/50 - Loss: -0.066   lengthscale: 0.202   noise: 0.030\n",
      "Iter 49/50 - Loss: -0.068   lengthscale: 0.196   noise: 0.031\n",
      "Iter 50/50 - Loss: -0.070   lengthscale: 0.190   noise: 0.032\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iter = 2 if smoke_test else 50\n",
    "\n",
    "\n",
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "for i in range(training_iter):\n",
    "    # Zero gradients from previous iteration\n",
    "    optimizer.zero_grad()\n",
    "    # Output from model\n",
    "    output = model(train_x)\n",
    "    # Calc loss and backprop gradients\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (\n",
    "        i + 1, training_iter, loss.item(),\n",
    "        model.covar_module.base_kernel.lengthscale.item(),\n",
    "        model.likelihood.noise.item()\n",
    "    ))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we reevaluate the model on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAADGCAYAAADSWvcAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7FElEQVR4nO3dd3iUZfbw8e89M5n0XkgISegJAULvKE2QRUBZLBR/vrZFYQUVWbfoioq6FlzLIijqWhFRUUREBUQWEVDpvUMoCamkkTblfv9IISGTSSBDSOB8rovrysw88zxnAjxn7nZupbVGCCGEMFzuAIQQQjQMkhCEEEIAkhCEEEKUkoQghBACkIQghBCilCQEIYQQgAsSglLKQyn1m1Jqu1Jqt1LqKVcEJoQQon6puq5DUEopwFtrnaeUcgPWAQ9qrTe6IkAhhBD1w1TXE+iSjJJX+tCt9I+sdhNCiEbGJWMISimjUmobkAqs1Fr/6orzCiGEqD91biEAaK1tQGelVADwlVKqg9Z6V8VjlFKTgEkA3t7e3eLi4lxxaSHEBdi8eXO61jr0cschGqY6jyFUOaFSTwD5WuvZ1R3TvXt3vWnTJpdeVwhRM6XUZq1198sdh2iYXDHLKLS0ZYBSyhMYCuyr63mFEELUL1d0GUUAHyiljJQkmM+01stccF4hhBD1yBWzjHYAXVwQixBCiMvIJYPKQoiry+bNm8NMJtM7QAek4kFjYQd2Wa3We7t165bq6ABJCEKIC2Yymd4JDw9vFxoaesZgMMi6o0bAbrertLS0+NOnT78DjHZ0jGR2IcTF6BAaGpojyaDxMBgMOjQ0NJuSVp3jY+oxHiHElcMgyaDxKf07q/a+LwlBCNEoHT582G3IkCGtYmJiOkRFRXW46667ogoLCxXA66+/HnzHHXdEX+4Yz+fl5eVwAo7RaOwWFxcX37p16/axsbHxM2fObGKz2Zyea//+/eY333wzyJXxSUIQQtSLxMREtx49esQeP368zmOXdrudm266qfXo0aOzEhMTdx09enTX2bNnDQ8++GCkK2J1xGKxXKpT4+7ubt+3b9+eQ4cO7V69evWBlStX+s+YMaOps/ccPHjQfdGiRZIQhBCNz2OPPRaxefNmn3/84x9Ob3S18c033/i6u7vbH3zwwQwAk8nEm2++eWLRokUhubm5BoBTp0659ezZMzYmJqbDI488EgGQk5NjGDhwYOvY2Nj4Nm3atH/77bcDAX7++WevHj16xLZv375d//792yQmJroB9OzZM/buu++O6tChQ7u//e1vEU2bNu1Y9s09JyfHEB4enlBUVKR2797tfs0117Rp3759u27dusVu3brVA2Dfvn3mzp07x7Vt2zZ+2rRptfrckZGR1nfeeefYe++9F2a329m/f7+5W7dusfHx8e3i4+PbrVy50rv09xm5adMmn7i4uPinnnoqrLrjLoTMMhJCXFLu7u5di4uLVdnjBQsWhC5YsCDUbDbroqKiLRdzzp07d3p26tQpv+JzQUFB9oiIiOI9e/a4A+zYscN7586du318fOxdunSJv/HGG7OPHDliDg8Pt6xZs+YQQEZGhrGoqEhNmzYt+ttvvz3UtGlT69tvvx04Y8aMyM8///wYQHFxsdq1a9degG3btnktX77cd9SoUbmLFi3yHzBgQLa7u7u+9957Y+bPn5/YsWPHotWrV3tPnjw5euPGjQemTJkSfe+996Y98MADGf/6179qXUMqPj6+2GazcerUKVPTpk2tP//88wEvLy+9c+dO9/Hjx7fctWvX3mefffbUyy+/3OSnn346BJCbm2twdNyF/F4lIQghLqkDBw7snDp1arOVK1cGFBYWGjw8POzDhg3L+s9//nPiUl63f//+OeHh4TaAG2644cyaNWt8brrppuzHHnssavLkyZE33nhj9vDhw/N+//13j4MHD3oOHjy4LZR0R4WGhpb3D40fPz6z7OdbbrnlzMKFCwNHjRqV+9lnnwVNmTIlLTs727B161afW265pVXZcWUJcMuWLT7ffffdYYD77rsvY9asWc0u9HMUFxere+65J2bPnj2eBoOBxMRE97oc54wkBCHEJRUTE2Px9fW1FRUVGUpbBQZfX19bdHS09WLP2aFDh4IlS5YEVnwuMzPTkJycbI6Pjy/69ddfvUr27jpHKUVCQkLRli1b9ixevNj/n//8Z+SqVatybr311qzWrVsXbNu2zWENNl9fX3vZz+PHj8+aNWtWZEpKinHXrl1eo0aNysnJyTH4+vpa9+3bt8fR+y9mNtaePXvMRqORyMhI64wZM5qGhYVZFi9efNRut+Pp6dnN0XueffbZJrU5zhkZQxBCXHJpaWluEydOTFuzZs3eiRMnpqWmprrV5XyjR4/OLSwsNMyZMycYwGq1MmXKlKhbbrklvewGvm7dOr+UlBRjXl6eWr58ecCAAQPyjh075ubr62ufMmVK5vTp009v27bNKyEhoTAzM9O0atUqb4CioiK1adMmD0fX9ff3tyckJJy97777oocMGZJtMpkICgqyN2vWrPi///1vIJS0MDZs2OAJ0LVr17y33347CODtt98Ors1nS0pKMv3pT3+Kueuuu1INBgPZ2dnGiIgIi9FoZO7cucFlYxj+/v62vLw8Y9n7qjvuQkgLQQhxya1YseJw2c99+vQ5XtfzGQwGlixZcmjSpEkxL730UoTdbmfw4MHZr7/++qmyYxISEs6OHj261enTp80333xzxrXXXpu/ePFiv7///e/NDAYDJpNJz507N9HDw0N/+umnh6dNmxadm5trtNlsavLkySndu3cvdHTtW2+99czdd9/dctmyZfvLnlu4cOGRP/3pTzEvvPBChNVqVWPGjMns06dPwdy5c4+PGzeu5auvvho+fPjwrOo+T1FRkSEuLi7earUqo9Gob7vttoyZM2emADz00EOpY8eObfXpp58GDx48ONvT09MO0LNnzwKj0ahjY2PjJ0yYkF7dcRfC5fsh1IbshyDE5eGq/RC2b99+rFOnTumuiEnUr+3bt4d06tSpuaPXpMtICCEEIAlBCCFEKUkIQgghAEkIQgghSklCEEIIAbggISilopRSPyml9iildiulHnRFYEIIIeqXK1oIVuARrXU80Bv4s1Iq3gXnFUKIaimlut14440tyh5bLBYCAwM7DRo0qPXljKsxq3NC0Fona623lP6cC+wFLlkJWiGEAPD09LTv37/fMy8vTwF89dVXfk2aNLl0NaqvAi4dQ1BKNQe6AL+68rxCCOHIddddl/35558HACxcuDBo7Nix5YXocnJyDLfcckvzjh07tmvXrl38xx9/HABUW0562bJlvj179owdPnx4yxYtWrQfPXp0C7v9ghf7NmouK12hlPIBFgMPaa1zHLw+CZgEEB3d4DYyEkJcpLvvJmrXLrxcec4OHcj/73+psRrq//3f/2XOnDkz4rbbbsvau3ev1z333JOxfv16H4B//OMfEYMGDcr5/PPPj6Wnpxu7d+/ebvTo0TnVlZMG2Lt3r+e2bduONG/e3NKtW7e4lStX+lx//fV5rvxsDZlLEoJSyo2SZLBAa/2lo2O01vOB+VBSusIV1xVCXN169epVcPLkSfe333476Lrrrsuu+NqaNWv8fvjhh4DXX389HEqK1h06dMgcExNjqa5MdMeOHc+2atXKAtC+ffv8w4cPm+v3E11edU4IqqTG7LvAXq31v+sekhCiManNN/lLafjw4VkzZ86MWrFixf7U1NTye5rWmi+++OJQp06diioeP3369GrLSbu7u5d/WTUajVit1so1tK9wrhhD6Af8HzBYKbWt9M8IF5xXCCFqNHny5PQZM2Yk9ezZs6Di84MGDcp5+eWXm5SNA/zyyy+e4Joy0VcqV8wyWqe1VlrrBK1159I/y10RnBBC1KRVq1aWxx9/PPX8559//vkkq9Wq4uLi4lu3bt3+8ccfj4SSctILFy4Mjo2Njd+3b5/HxZSJvlJJ+WshriJS/lpI+WshhBA1koQghBACkIQghBCilCQEIYQQgCQEIYQQpSQhCCGEACQhXFbJyckMGDCA06dPX+5QhGiUjh8/bho5cmTLqKioDu3bt283YMCA1jt27HCv+Z2Vff/99z6tW7duHxcXF3/06FG34cOHt3R0XM+ePWPXrl3r0rpNDYnLituJCzdr1izWrVvH008/zdy5cy93OEJctFnL9jR15fn+OTI+qaZj7HY7o0ePbj1hwoSMZcuWHQHYsGGDZ1JSkltCQkJRTe+v6MMPPwyaPn168pQpUzIBvv/++yMXF3njJi2EemK12UnPKyIpqwAPD0+UUsybNw+73c68efNQSuHp6Xm5wxSi0Vi2bJmvyWTSjz76aFrZc3369CkYNmxY3n333desTZs27du2bRv/9ttvB5Yd76i89b///e+Qb7/9NujZZ5+NHD16dIv9+/eb27Rp0x4gLy9PjRw5smXLli3bDx06tFVhYWF5baMvv/zSr3PnznHx8fHt/vCHP7TMzs42AERGRnZ8+OGHm8bHx7dr27Zt/NatWz0AsrOzDTfffHPztm3bxrdt2zb+/fffD3B2nstBEsIloLUmObuA9YfS+WrrSd5dd5Q5Px3iow2JLPr9BH9/fyVdB43Ezd0DADd3D64ZfhMLVv3O0fSz2OxSDFaImuzYscOzU6dO+ec//+GHHwbs3LnTc+/evbt//PHHA0888USzxMRENygpb/3GG2+cOHTo0O7jx4+7r1y50mf69Onp1113XdYzzzxzcunSpUcrnmv27Nlhnp6e9iNHjux+5plnkvbs2eMNkJycbHruueci1q5de2DPnj17u3btmj9r1qwmZe8LCQmx7tmzZ+/dd9+d9vzzzzcB+Nvf/hbh5+dnO3DgwJ4DBw7sueGGG3JrOk99ky4jFzqekc/+lFyOpudxtqj6gll+wWF4ePlgLS7CZHbHWlyEzeRJYr6ZxK2n8HAz0jrMh7ZNfIgO8qKkoKwQojZ+/vln31tvvTXTZDIRFRVl7dWrV966deu8/P397Rda3nrdunU+06ZNS4WSUttt27bNB1izZo334cOHPXr27BkHYLFYVLdu3cr3TZgwYcIZgJ49e+YvXbo0EGDt2rV+n376aXlXVGhoqG3hwoX+zs5T3yQh1JHFZmdfci7bTpwhPa/Y4TE5Gal8+Nx0xkx5jK/mPssdj71CblYGfUeOp/eI29i4fBE5meWtXgotNnadymbXqWyCvM30aB5EXLgvBoMkBiHKdOzYsWDJkiWBF/IeV5W31lrTv3//nG+++eaoo9c9PDw0gMlk0s6uUdN56pt0GV2kIquNDYczeHfdUVbtTSE9r5icjFTmPHJ7pZs7wIoFczm6axMfPz+Do7s2seLjN7hr5hzGTp1JZKs4xk6dyV0z5zi8TubZYn7YfZoPNhxj16ls7NKdJAQAo0aNyi0uLlazZ88OKXvu119/9QwICLB+8cUXQVarlaSkJNNvv/3mc80115y9mGv0798/b8GCBUEAv//+u8eBAwe8AAYOHHh206ZNPrt27XKHku06a5rdNGDAgJxXXnklrOxxWlqa8WLOcylJQrhANrtm6/EzvP/LMTYeyaCg+FzXUNmNf8XHbwDw6MgEpg+LZf2yhWitSUk8hNaa9csWMn1YLI+OTKjxemVJ5vjJJFbuSWHBb8dJzi6o8X1CXOkMBgNLly49vHr1ar+oqKgOrVu3bv/Xv/418s4778xs3759Qbt27doPHDiw7VNPPXUyOjraejHXmDFjRurZs2eNLVu2bP/YY49FxsfHnwVo2rSp9a233jo2bty4lm3bto3v3r173M6dOz2cnetf//pXclZWlrFNmzbtY2Nj45cvX+57Mee5lKT89QU4mJLLukPpZOVbKj3/6MgErMVVZ7kZTG50vuZ6dvyystLrJrM7Cf2HMXrSX/ELCnV6zS9ef5IN335KnxvGcfO0JwFQChKa+dO3VQgebsY6fy5x9ZDy18JZ+WsZQ6iF7HwLP+1P5Wi641bn4x+sYun8F9i5fhWWokKUwYi22+g66AbMZg9slmKUwYC221EGAzZLMR5ePk6TwaM3dMRqOTcmsX7ZQtYvW4jJ7M6Ly3aw/UQ2h1PPMqRdGC1DfVz+mYUQVx9JCE7Y7JpNxzL5/VgmFlv1LamyWUOWokIAtL2kG2nTyiUAKIOBVgk9yTuTgU9gME2iWlYZZzhf5wEj2LRqCQaDEbvdhpu7Bx37DWXQLfcw55HbueOxVyAolK+3JdElOoBr2oRilEFnIUQdSEKoxunsQlbuOV3tzKHz5WZl0GPoGHKzMjiw5RfsNhsmszsent7c98J7RLaMq9V5zu9+spcmF0tRIR5ePmz4dlH5OEVZF9LW41mcyirgho4RBHg5nUUnhBDVcskYglLqv8BIIFVr3aGm4xvyGILVZmfDkQy2JGZhP+9342j66PndPl+8NpMNyxdhMLljs7QDOtG01RiCwgaTdsqN4kIDSoEyagwKfIOsRLQoJjAsjS2rn+fmB29i3ZL5lbqfYrv1Y//mdWh71a1fy7qQAMwmA0Pjm9C2ie8l+/2Ixs2FYwhHOnbseMZgMMi0t0bEbrernTt3Bnbq1MlhrSZXtRDeB+YAH7rofJdFUlYBK/ekkHm2cqugLBEEhUeWTx9NPX640rd0gOIixakjrdD6TWyWkUBJeZekw0UkHT6EUgfpPnQQ2q6w28Fug6w0Nzat8qUoPwB4lznTCwlqEoylyITRbTV2ay5BTSKZueB/lcYpyrqQRk/667nrW+18uyOZ9BZF9GkVLAvaxKW0Ky0tLT40NDRbkkLjYLfbVVpamj+wq7pjXJIQtNZrlVLNXXGuy8Fqs7P+cAZbjp/BUYPpqYkD0HY7R3b+DkBK4iHg3ECv0S2O/qPW8+v3fhTmP4PZw4qX70ZyzzyOzboWk/k0Cf0Hl84qSql07nNdRDFAJ+y2EaQnjQW+wmi0ENJ0LRnJHzhc3ezh5QNal48p+AWFkpORyv89cjtPvDKfcQM64m6SWUjC9axW672nT59+5/Tp0x2Q6euNhR3YZbVa763uAJdNOy1NCMuq6zJSSk0CJgFER0d3S0xMdMl16yo5u4AVuyu3CspaBMf3ba800+d8RtNgfAKeITujNwYDdLomlx7DcmidUMCSeU+wYfkijG5mbJbiStNGK8rJSHUwQ0nRpsvjBIQ8xJY1vmibovvQHLLSHiQ00lJpdbNvYEilaakVp6ne/4/nGN0pEn8vt0vwmxONkau6jMSVqd4GlbXW84H5UDKGUF/XrY7VZmfjkUw2J56pMlZQtsCs25Absdus5TfrcxKAF7FZr+dsdi5Dbsuk36hsAkLPrX1xVpqioupmKB3c+iTwJEa3FvQbvZkNy/2xWd8jIDQH/5A0fv3+C4fTUs9/bDK7cyI1i3D/y7bWRQjRSNRbC6Gi+h5UTk5OZty4cSxatIjw8HBOZRWwysFYQXULzACUMqB1M9zML2EpvhmDMY/o2G/w9vuUe55+uU7xvffUA3h4+VSaoVRxjMAvKJScTCM/fRbEz18H4OFlJyxqLsf2TCuflmoyu+MbGELumXSsxUWV3h8cGsbwDhG0DpP1Clc7aSEIZ66Kvr+yjWhmPvkUP+1P5fNNJ6okAyhZYFaxLLUylPS/h0S2pc/I3SjDQayWmxh0axZPf5bKtFd71DkZANw1cw7j//I8QWFN0XZ7pTGC8llMOpkTB69HqS7k527g2J6pwE/Y7a0BsBYXYXb3xGYprvJ+i02zbEcSW4+fqXOsQogrl0u6jJRSC4GBQIhS6iQwU2v9rivOXReenp4UFp7r6pn/1pvMf+vNSlM1K3LcfTOQ9FPzSD8VhzJ8zuMfdSMw7KLKotTIWTdTxW4sm/VFdqxbgM36LLCVJjFv0KzNeg5uWVft+7WGNfvTyC20ck2bEJmBJISo4oquZZScnMy0h6azdOkSigsLq3TDOFLWfZOVbuHgtltA3wEcoXWnD7j978NqrD3kas66sYxuUdgs84HhdLoml1sfTsHTp+pahfO1i/BlaHy4rGy+CkmXkXDmiu0ystjsJBaYSS5QWIqKHHfDOHDXzDm07zOHY7s/Bj0eZXge6EBY1MEq76uu3LUrnd+N5ebugX9ION2HjuGh1+fQd+THhMe8zc71PsyeHMOxvTUPHu9NzmXp9lMUW2tOHkKIq8cVV7rCZtfsOpXNr0czOFtkIyczvUo3Stm00vNXGicfzWTuozmczR6Gh/dB2vd5lyHj+rJx+RiHN/2ybpxl78wmM+WUw5XLFSkFRqXQlHThaLTDdQ8VOVp/0L73oPIprDdPewKAxL0n+Oi5CN54JIo/PpBKnxHZTs97LD2fxVtOclPnSDzNslZBCHEFdRnZ7ZoDqbl89+te5vxzqtObs6OS0tt/9mHBC95Yiz2Jjv2Kqa90wVhNuqyuG0cZDLz8/V6gpIxEdJAXTfw8CPJ2I8DLTICnGyZj5UZZXpGV9Nwi0vOKSM8r5uSZfHILK49RvPfUA/gFhVZKao421CnIM/DRcxHs2+RNv9FZ3HR/arWfoUyQt5mbukTi7ylrFa4G0mUknGn0CaHQYmPnqWy2n8git9Dq8GZfxvGNPAj4DzAB2ATcCewGqHbwuWwx2ZafljmMyd3dnbyz+VVu/rWhtSYpu5ADp3M5mJrrdG9mR+w2WPZuCGu+CMLDexNTX7ER0dz5LoO+HiZu6hJJiM9l26hJ1BNJCMKZRpkQbHbNqTMF7E/JZf/pHCw2Xe239oo39fNXBaNuBD0PpUIIb/4paaemYC3Oq9Xg8xevzWT9t5+Wrk8o6Yv39PLij2PGMHv2bMLDwy/685XRWnMoNY9fj2aSlut4YLk6cx9dxaFtd2P2zCKs2VTunfWQ0+4sDzcjozs3JTLAs65hiwZMEoJwplEMKlttdtLzitiTlMO3O5J583+HWbzlJLtOZZfvU+Bo8LXr4FE8/uGPlc51cPuvWIq8gA9BLwFS0bo7yUfvwFqcV+vB59ysDIbfegcjRt8IlGznV1RYiJ+fX43JIDk5mQEDBnD69GmnxymlaNPEl4m9ohnVqSlN/KofMC4b4H70ho5MHxbLoW1/BgZQXKA4efAtnhz/WKXjzh8TKbTY+GrLSY6k5TmNSQhx5WrQCeGb7Um88/MR5vx0iI82JPLD7tMcSMl1ODumuuJvFW/qKxbMJTfzWpRhP6jxoGYBPVCGki6i0GbNefC1z+g7cjy5Z6rfHTDC34Pl3yzhu0UfYDbAlClT2LJlC/fff3+NN3k4t1Du6aefrtXvQSlF6zAfJvSK5g8dw/F2rzoIXDbA3XnACLoOGln67G9AbyAV9A9MH/YET00cUGnf54oJwmLTfLM9mV2nnA9ICyGuTA26y2j+2sMX1Ide3eBrSXdSEPA6cDMlYwX3AFXHB6D6sYMALzf6tw6hzUXuN3D+QrkyHh4eFBQU1Po8hRYbvxxKZ+epbP5yQ03lNuxAILAEuBZ4FHipwkEGFLrKmEvfVsH0ahlc65hE4yBdRsKZBpsQkpOTGfCHG5n4d+dTOWtis8GqT0ys/CQcu80AzMRkfgNvP19i4rtQXHC22vpBZUwGRffmQfRoHnhRA8UVP9OMGTNYsmQJ+fn5eHl5MaYOYw5JWQV8vnYnH746q9I+CV6+AbTp0ofiwgJ2/Px9aWJww+z5GcUFo4HXgIeBqn/3FZNhpyh/BsWGyarmK4gkBOFMg+0ymjVrFod2nOvaKHMhi8GO7vbg33+O4YePWuIXtB/ogMn8GjZLPu17D+LOx191Xj8IaBboycTeMfRpFVynZAAQERGBn58fhYWFeHh4UFjLMYfqNA3wZPKI7kSEBlVZpzDhL8+j7Tb6jZrA9Llf0W/UWLz9/gz8G3gQWABU3m7z/DGX7SeyWbYjGYtNFrAJcTVocAvTzu9WqVjG+cVlO8r7ys/frayi9CQ3vns/mK1r/AgItXDnE0lsWvVX2vfuRe8R/65U56e6+kHubgaubRNKh0h/l36+lJQU7r//fiZNmsT8+fNJTk6u0/nMJgNmSy7j7ribmP438vPST8s/Q8W1CmOnziTnqQdo1XEdifs8STs5GQgB/ogy5KPtdg5t21jl/IdS81i8+SSjOzfFy9zg/rkIIVyowXUZnd+tUtaNs+PnHxxuVmN0MxMT14k7HnsFpcJZsSCYDd/6YzJprv3jGYaMy8Td88I+Y8tQb4a0a4KPe+O6AWblF7NsR7LDKaoVV2ev+GgO6781A+8C24nvPZvM5E2cTjxE35HjHSbaAC83buocSaC3ucprovGQLiPhTIO741XsVnGr0I3z+Ic/OtxT2GA0smnlZub/4wzpSX2wFit6j8hm2MQM/IIvbFGXu5uBgW3DiG/qd4k+3aUV4GXmth5RrNidwoGU3EqvVWxZ5WZl0G9UKGHRm/jqjQT2bHwaGAZUbZGVycq3sGjTCUZ1krUKQlypGlxCgHPdKkHd/sBPSxaSk5lWZVqppaiQLatPAI8An5F0RAGfYnR7jpunfXnB12wR4s118Y2vVXA+N6OBGxIiCDliZsORjCqzkMp2VSu54c8kKOwg789qic26HhiOm/sBOvYbyqBb7qm0VzNAQbGNLzefZEi7Jo02aQohqtcgB5W//PJL3njjDaJat2PohMmczckiJzON3KwMug+9j2vH/IaX715gK3Aj8B9M5vZ0Hfwp//zorQu6ltlkYGh8E27qEtnok0FFvVoGM6pTU55asNrpgr32fTzp2PcpwAqsxVLUg0PbNvK/L9+vtF6hjNWu+WH3adYdTOdydDcKIS6dBn8H/OHjuRzZmcF7T20n9cRMCvP7ou0KL99ESloH72IyF2KzFOPh1e2CpqhGB3kxtH0T/DyuzMJurUJ9uPf6rqz42M/pgj2bbQfdr3uRI7ueI/P0CnIyx7Np5RKg+i6k349lkplfzPD24ZhNDfJ7xVXjWPpZmod4X+4wxBWgwSYEs3ksFssNwHNAcxL3AhwHXgAWkJ+7q/xYa3FJpVFnq4srndtUMoOoYzPXziBqiMJ8PfDjLIPHTKTz0Fuq7KQGVFi89wmwDPgC+DNQ0trqOngUoyf9tcq5D6fmsajgBKMTmuLvdWUm1YYsO9/CmgOpHEk7y0PXtZH1IqLOGmxCGDPmIxZ/acNmXUHJytofgf1VjqtNIbqKWoZ6MzguDN8rtFXgyNIlSyiy2li2PZnIVnEOj3n8g1Usnf8CO365AWvxB8CbQCQw02ldp/TcIj757TjDO4TTQr6l1guLzc7vRzPZnHgGq1267YTruKStr5QarpTar5Q6pJT6myvOOX++F33/cD9wM0q9SVkycHP3IKRpDEqpSl0gaO10wZq3u5HhHcK5sXPkVZUMyribjNzUJZJ2EY4Hg8sG7W2WLGAsJVNS/0mT6B/Jycx0eu5Ci42vt51i45EMGVe4xA6l5vLB+mP8ejRTkoFwuTonBKWUEXgD+AMQD4xXSsXX9bz+/pCblU6/URPo2H9Y6bUMWIuLsNtt9B05vlIhuorTKivHB52jArijT/Nqb4ZXC6NBMbxDON2r2R+hbJHeI/O+oO/I1YQ2W0DK8UFYLYspzD/XHeFotbjWsOFwBku3J1FoubDpvqJmWfnFfLX1JN9sT66ygZIQrlLnhWlKqT7Ak1rr60sf/x1Aa/2v6t5zocXtnO0Y5mwfhI9+PsDguDDCnJSNvlptOX6GtQfSatzCc8O3/iz+TxgRLYu4d9Yp/INtVTYhOn9LUl8PEyM6RtBU1ivUmdVm57djmWw+5rx7qLZjCLIwTTjjioRwMzBca31v6eP/A3pprR8477hJwCSA6OjobomJiTWeuzbVTs/f9MbN3YNu117PS7Nn069jKxloc2Lf6RxW7E7BVkPXw97fvPjgmaYUF56ipCG4u9LrymAAXbliqkEp+rYOpntMoPwdXKTjGfms3pfCmXxLjcdKQhCuUG/zBbXW87XW3bXW3UNDL7566fkqLlgrWdlcTKdWEfRPaC03ohrEhftxY+emTqeN5mSk8uOiP3LXzB14+4cA6zGYRlQ6RtvtaK1Zv2wh04fF8ujIBOxas+5gOku2neJskXRxXIj8Yivf70pm8ZaTtUoGQriKKxLCKSCqwuNmpc/Vm/ycTMZOvItf1q9n8uT7SU1Jqc/LN2oxwd6M7doMT3PVTXfgXMmLnb+8SFz3R4Fj2K1fYzA+CEBI05jyRW8mszs+/kE8+Ppn5e8/lp7PRxsTOZSa6+j04jy7k7L5cEMie5Nr//vKyUhl4MCBtdqcSQhnXNFlZAIOAEMoSQS/AxO01rure4+rNsjx93Sje/NA4iP8qpSmTk5OZty4cSxatMgl+xtf6TLPFvPV1lPkFJR8I61ubAZ8gI8pWSH+JoFNnicr9ThGN3P58dUVyItv6sfA2FDcTY6Tz9Usu8DCj3tTSMzIv+D3fvH6k2xcvoj77ruPuXPnOj1WuoyEMy6pdqqUGgG8ChiB/2qtn3V2fF0SglIQE+xFQrMAWgR7YzA47haaMmUKb731Vq3+k4gSuYUWlmw9RXpecZWxGWUwAAptt2EyexEc/iEpx8diMP4Pu/1m0FUXBRrdzLz07c5Kz/l5ujG0XROig73q6VM1bFprtp7IYsPhDIdbw1ZUNng/ZspjfDX3WRL3bsNmtQB9KNkVbzngfAc+SQjCGZeMIWitl2ut22qtW9WUDC6Wl9lIt5hA7uzbnDFdmtEq1MdhMvD09EQpxbx587Db7cybNw+lFJ6eMuOlJr4ebtzSPYrIQM8qxQS13V6aDNyxWQpolbCY1p3mYbf1xuy+i7juU8u7jpShpAXQZeCIKtfIKbCweMtJVuw+fdVPT03PK2LR7yf4Zv1u/v3ghBo3fSrrvvv4+Rkc2bkJb//xwBpgPfA0np5eTJw4kaNHj9ZD9OJK1KCL0BiUolWYD6M6NeXea1pybdtQAryc1+M/cuQIEyZMwMur5Buol5f8J7kQHm5G/tglktZhPuXrEh587TOCwpsRFN6MB1/7DJRi/bKFHNo+BbiW4kIL+zY9j6VoLADaXnKj37RySfkg8/l2J+Xw0Yarc2zBZtdsOJzBJ78eJzm7sNo1NGUeHZnA9GGxrF+2EK01KYkdgU3kZHwAtAQeAq6loCC/TjvwCdHgNsipqNBiw8PtwvubJ0+ezPz58zGbzRQXF0u30UXQWrN6Xyo7TmZXee387iSTORptX4jN2pfAsC/ISr8dbS+qdVmRFiHeDIytOdlfCZKyCvhxbwrpecXVjtNU3PTJLyi0/Pe9/ed8bNZnKOki2k9JXa+PUQYbw6+/nrCwMHJycvjyy+rLv0uXkXCmQbcQLiYZwLn9FDZu3Mj9998vsy8uglKKIe2a0LdVcJXXzu9OsllO0GPYPAaMzeRM6s1o+08Y3do4rKzqyNH0s3y0IZH1h9Ov2P2bi6w2Vu9L4bNNJziSeJI5j9zOPU/Pw9s/CJPZHThXmrzLwBGVWgwFZyM5tucJbNYfgWjgbkqKAryHMthAa2JiYnj//fedJgMhatKgWwiiYdiTlMOqvZUXsFW3evzVqYs4dfgRTG4mWrSfg5v795X2dgaqrGyuyM/Tjf6tQ2jbxKdRryOpOMstV3mzZn8aeaXrMcpWeodFtyIl8RBAeV2uynyAJ4CHMBiLiWqzDKN5HvnZSeRmpeMbEIJPYDBNoloSZiqoVTKQFoJwRhKCqJUTmfl8syOJIkvN3+Azkt348NkIThzwwC94EdNejSKoSUj56+eXvnAk3N+D/q1DiApqnLORyma5XffHiQyf9DjgbCrvOUaTGygDNsuNwCtAU7oMSmHM5Dx8AqofhJeVysIVJCGIWsvIK2LJtqTytQrOWIsVr0w9SPLREXj5JjJltubVqXHV1p2quPlORS1DvenTMrjR1KPy9PSksLCwyvNGNzP/LN0XfMcvKyv9HkxmdxL6DyvdH3wnJaXHrwO20LHfEsY+MKzaFlUZSQjCFRr0GIJoWIJ93BnXI4oIf+c350dHJvDoyLYkH70BGEl+rgez74vAZp1Kl4Gjqt3O05EjaWdZ8Otxvt52iuRsx3PrGwqrzc6StZvpMWSUwym450qMF5eu6yipA2UtLmLL6u/YtDIW2An0AKaA6gVsrHEWkhCuIi0EccGsNjur9qZUW17h/FlIqCag3wRuwi94NzkZN2Ayn8ZmKXbabeRITLAX3WICiQ7yKv9GfLlXpdvtml1J2fx2NJPcQitfvDaT9d9+6vBYpQz0HTmOlBNHyDuTgU9gMG7ma9i/+Q7stg7AYkzmv5DQvwM7fv4Bq6W4yjkctaikhSBcocHumCYaLpPRwPAOEQR5u7P+cHqVEtpl34QtRaVdJzoFGAPcRU7Ga8Aurh1zhMKzL5J7JrXKClxnXSOJGfkkZuQT6OVGQlQA8RF+zJo1i3Xr1vH000/X6/TiQouN3Uk5bD+RRXaFbrTcrAx6DB1DblYGB7b8gt1mczgFt7hQ8f2Hwfzvy0Dc3M5QbBuDyfxd6f7gvXm8tIupYiXfsnMIcSlIQhAXrWeLIIK8zfyw+3SVsguOb4oLadtVUZg/m9WLEoiKfQNtv5dl/325fAVu6vHDrPj4jRpbDWfyLfRqE1GpL37evHnMmzfPaekGV8g8W8z2E1nsSc5xWG6ibFbVF6/NRNvtlXb2K0sGB7d68tmrTchINtPnhiyy0+8mMMyb3iM+K5+1df703tpO4xXiYkmXkaizzLPFfLM9icyzVbs3vnhtJhuWLyovfucXFMrDb3zFgS0tWfiSO+ANPAe8CFQejD2/a+T86aqO9sK45roRPPf8i3SOa+7SInrpeUUcSs3jYGoe6bnOZwqVcTQ197bpc/nmnVB+/c6fkKbF3Do9hdYJ1ScvZ5tDAQR5m4kN96V3y6rrRRyRLiPhjCQE4RJFVhs/7E7hcGpepecr3tAWPP8IpxMPoQwGtN0OhAGvAeOAw8CDwLfls24G3XIPX819trwrKSg8ks2rvq407lAx4VQckzAoRaivO80CPYnw9yDAy0yAlxtuxprnUVhtdtLzijmdU0hKTiHJWQV13pdAa9j2Px+WzAvjbLaRAWPPcP0dGZjdL/z/n4ebkXYRvsSF+xFewwD/+SQhCGckIQiX0Vrz+7EzbDicgb3Cv6ua598PBv5Dyerbb4CH6TuyJwAbvv0UXXLyKu8ymd1p1+Nap9+gK1IKfNxN+HqYMBkMmIwKk8GAUlBQbCPfYqOg2EpBsb1S/I44W1x3vswUE4v/E8be33xo1qaQWx9OoVnr2rUyKvIyG+kaE0hCM/+Lbv1IQhDOSEIQLnfyTD7f7zpdvhm8o64d/+AmZCQfBxRa2/ELikTbp5KXMxVtNwFzgWeADIfX6Dp4VI01ki6l2iyusxYr1iwOYNUnwWg0fkGvMWV2LwJruWNgWdK578nXGdItloRm/rVq4TgjCUE4I+sQhMs1C/RiYq8YWoR4A1VrH1mLi7DbbfQdOZ7pc7+i36gJxLSL56nPxvDEgmS6Ds4EplLSjfRX4Fy3iDIYUEpd9OBqTkYqcx65vcZS09U5v/JoxW1DK9q90ZsXJsWw/L1QYrufpWO/h8k8PYMfF9Z+LcHKT+ZydNdmDq94n24xgXVOBkLURFoI4pLRWrPl+Bl+OZTBOzP/XOuuHYAPn5nPtrXXAqOBZGA28Cadrr0WH/8gcjLTGPvAE9V221TXpVObb/bOOGrtVJxOenBrHu8/k09Bbn/CoopIT/4jduvyKudxtjq7ui42V8yekhaCcEamnYpLRilFt5ggooK8CHxpPul5JbOQxk6dWeN7bbYd9BuVzLG935F28g6KC1/GZJ7J6cQvePDVa/HwtvPF60+Wr+A9/+ZecXXvzdOerHKTXb9sIeuXLXR6Y3akuqmgRQVN+ehfwWz9yQfIoXn8x/x5dk/OZv+dpfMNtV5L4OthYu3mXcz510yWLFlCfn4+Xl5ejBkzhtmzZ9c6TiEuhiQEccmF+Xowvmc06w6ls+1ElqPx4SrObz0c3X2clZ8Es+/3u/nHmBzgQ2AroCvd3NG60uresteMbma6DhrpkkVeZRsH9R5xG6sXrWXvb0NYvywKKAKeB2ZzbM8Z/jKipCXQc+iYGtcSKAWdogLo2yoYd5ORD/38KCwsxMPDg8LCQtn4RtSLOnVKKqVuUUrtVkrZlVLSDBXVMhkNDIwNY2zXZvh7ul3w+1u0L2TSs6d4eE4iwRGbgT8Be4FVGE0T6TxgLI9/+COdB5Rs22korSFUVi/podc/4+D2X8tvzJaiQg5t23hRn+X//XMO7Xq9yLJ3h7B1zWxyz1xPr+GpdOx3P27us4Azleo0Vdx5ru/I8eSeqbz/dIivO7f1iGJQbFj57CHZ00NcDnUaQ1BKtQPswFvADK11rQYGZAzh6lZstbP+cO1bC2Uqd/uEAPcCkynZNCYP+BpYBPwAnGsl9B05HihpLTSJac3tf3u5fE1E35HjuXnakzVOI9Uako6Y2bHOl60/+ZKeZMYv2Eq/UVn0/kM2voG2atdEVMdoUPRsEUSP5kEYHewPfinIGIJwxiWDykqpNUhCEBcoObuAlXtSyMirusLZkfMHdJXBSNuu16D1NaQc70HemQHYrH7AWeAXQiOPkZb0AejfAKvTcyuDAbSudBPPyTRy4oAHR3Z6suMXHzKSzCiDpmXHAvqOyCbhmlyMpnOxvXjfKOJ7DWLAH++sceC8aYAH17VrQrCPe+1+WS4iCUE4U29jCEqpScAkgOjo6Pq6rGjAIvw9mdgrhs2JZ/jtaAYWm/MvJ1W37iwmODyCm6fdDsBnr0xh43eZKMMNaPs1pJ0aBkxCKQuaw6D3owyH8fDOpCg/DbstFygANNoeAoSxflkY65ftpqQEdbPSK1tplZDL4FvOENPuGF++MZVWnV7BaDrXklixYC4FudmY3T2IbBVX7cC52WSgb6tgOkcFNOod4cSVqcaEoJRaBTgazXpMa/11bS+ktZ4PzIeSFkKtIxRXtLJuk7gIX9YeSONgSp7T4ysO6JZ9Cy9zNieVfqNC6T0iio3LXyAzpZie17/Gqk9+59RhC6hYtH0YBbmeTq6Qg7tXDm06u9OyQypHds5j14bZBIZdx+bVpzi6J7Ly7KUbOjocxHY0e6lVmA+DYkPx9bjwMRQh6oN0GYkGJTHjLEvX7+bVxx6oVVmI2qhYT2nDt5+xY91m8rILQLsDnpTMrUgr/VNSgC8/N9vhXgRVKAVaYzAYsdsdl7n2cTcxKC6U1mG+df4sdSVdRsIZWfooGpSYYG/2LH+Po7s289PCeXU+X05GKmdzshg6cQqRreK4edoTPP3Z1zz5ycd0HdQWN/f9wFaUIZm4Hr0IbdaCnMw04nsPxts/qOYLlH6hsttL9ju2FBWWTys1GhQ9mgfx//o2bxDJQIia1GkMQSk1hpKqZKHAt0qpbVrr610SmbjqnL8f8f++XsD/vl5wwYvHKjp/gVoZRwvM9v2+tvz1HT9/X/6zUga0PrfvQVm11pCmMWRnpJQPcMd264dPQDC5Z9JpGerNgLahBHiZLypuIS6HOiUErfVXwFcuikVc5Y4cOcKMGTOqrNB99vkXSLF4suX4GQqKbbU6V21WJlccj1j75fsc2LqB7PSq8/0rJgOAhP7D8PEPYu/v/6s0wB3UJJI/P/48fVsF07y0jpMQjYl0GYkGIyIiAj8HK3RjmkXSs0UQd/drwbVtQ/D1qPl7zOMfrKLroJHlm91XXChW5q6Zcxg7dSaRreIY/5fnad9rIFDSIijj5u6Bf0g43YeO4ZF5X9Nv1ATsNlvp++LLF5wNHjMRb3seE3pFSzIQjZaUrhANStkK3UmTJjF//nySk5PLXzObDHSLCaJLVCCH0/LYdiKLk2ccF3u7mO0nc7My6DdqArlZmez4+XuUMmAtLqJ970Hl3U0Vp5PeNXMOEf4edIkOpO19o2QaqWj0pNqpaNTS84rYnZTDwZTc8v0XytS0/WR1anqf0aBo28SXzlEBF7xj2eUms4yEM5IQxBVBa82prAIOpORyJO1sleRQV0aDIjrIi9ZhPrQK9cHT7Lr9muuTJAThjHQZiSuCUopmgV40C/RicBxk5RdzIrOAk2fyOZ1TSE6BtcZtMSsyGRRhfu6E+XnQ1N+TmGAvPNwaZxIQorYkIYgrUoCXmQAvMx2b+QNgs2uy8os5k28ht9CCza6x2DQ2u0aj8XQz4mk24mU24eNuItjbjKGeCs4J0VBIQhBXBaNBEezjXu/F5IRoTGTaqRBCCEASghBCiFKSEIQQQgCSEIQQQpSShCCEEAKQhCCEEKKUJAQhhBCAJAQhhBClJCEIIYQAJCEIIYQoJQlBCCEEUMeEoJR6SSm1Tym1Qyn1lVIqwEVxCSGEqGd1bSGsBDporROAA8Df6x6SEEKIy6FOCUFrvUJrXbYTyUagWd1DEkIIcTm4cgzhbuA7F55PCCFEPapxPwSl1Cog3MFLj2mtvy495jHACixwcp5JwCSA6OjoiwpWCCHEpVNjQtBaX+fsdaXUncBIYIh2skGz1no+MB9K9lS+sDCFEEJcanXaMU0pNRx4FBigtc53TUhCCCEuh7qOIcwBfIGVSqltSqk3XRCTEEKIy6BOLQStdWtXBSKEEOLykpXKQgghAEkIQgghSklCEEIIAUhCEEIIUUoSghBCCEASghBCiFKSEIQQQgCSEIQQQpSShCCEEAKQhCCEEKKUJAQhhBCAJAQhhBClJCEIIYQAJCEIIYQoJQlBCCEEIAlBCCFEKUkIQgghAEkIQgghStUpISilZimldpTup7xCKdXUVYEJIYSoX3VtIbyktU7QWncGlgFP1D0kIYQQl0OdEoLWOqfCQ29A1y0cIYQQl4upridQSj0L3AFkA4PqHJEQQojLQmnt/Eu9UmoVEO7gpce01l9XOO7vgIfWemY155kETCp9GAvsr0V8IUB6LY67nBp6jA09Pmj4MTb0+KD2McZorUMvdTCicaoxIdT6REpFA8u11h1ccsKSc27SWnd31fkuhYYeY0OPDxp+jA09PmgcMYqGr66zjNpUeHgjsK9u4QghhLhc6jqG8LxSKhawA4nA/XUPSQghxOVQp4SgtR7rqkCqMf8Sn98VGnqMDT0+aPgxNvT4oHHEKBo4l40hCCGEaNykdIUQQgiggSQEpdRwpdR+pdQhpdTfHLzurpRaVPr6r0qp5g0svulKqT2lZTx+VErF1Gd8tYmxwnFjlVJaKVXvM1JqE6NS6tbS3+VupdQnDSk+pVS0UuonpdTW0r/rEfUc33+VUqlKqV3VvK6UUq+Xxr9DKdW1PuMTVwCt9WX9AxiBw0BLwAxsB+LPO2YK8Gbpz+OARQ0svkGAV+nPk+szvtrGWHqcL7AW2Ah0b2gxAm2ArUBg6eOwBhbffGBy6c/xwLF6/h1eC3QFdlXz+gjgO0ABvYFf6zM++dP4/zSEFkJP4JDW+ojWuhj4lJIprBXdCHxQ+vMXwBCllGoo8Wmtf9Ja55c+3Ag0q6fYah1jqVnAC0BhfQZXqjYx/gl4Q2t9BkBrndrA4tOAX+nP/kBSPcaH1notkOnkkBuBD3WJjUCAUiqifqITV4KGkBAigRMVHp8sfc7hMVprKyVlMoLrJbraxVfRPZR8S6tPNcZY2n0QpbX+tj4Dq6A2v8e2QFul1C9KqY1KqeH1Fl3t4nsSuF0pdRJYDkytn9Bq7UL/rQpRSZ1rGYlzlFK3A92BAZc7loqUUgbg38CdlzmUmpgo6TYaSEkra61SqqPWOutyBlXBeOB9rfXLSqk+wEdKqQ5aa/vlDkwIV2gILYRTQFSFx81Kn3N4jFLKRElzPaNeoqtdfCilrgMeA0ZrrYvqKbYyNcXoC3QA1iiljlHSv7y0ngeWa/N7PAks1VpbtNZHgQOUJIiGEt89wGcAWusNgAclNYQailr9WxWiOg0hIfwOtFFKtVBKmSkZNF563jFLgf9X+vPNwGqtdX0toKgxPqVUF+AtSpJBffZ71ypGrXW21jpEa91ca92cknGO0VrrTQ0lxlJLKGkdoJQKoaQL6UgDiu84MKQ0vnaUJIS0eoqvNpYCd5TONuoNZGutky93UKLxuOxdRlprq1LqAeAHSmZ6/FdrvVsp9TSwSWu9FHiXkub5IUoG1cY1sPheAnyAz0vHuo9rrUc3sBgvq1rG+AMwTCm1B7ABf9Fa10tLsJbxPQK8rZR6mJIB5jvr8YsJSqmFlCTMkNJxjJmAW2n8b1IyrjECOATkA3fVV2ziyiArlYUQQgANo8tICCFEAyAJQQghBCAJQQghRClJCEIIIQBJCEIIIUpJQhBCCAFIQhBCCFFKEoIQQggA/j8l52LEGUwMGQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    trained_pred_dist = likelihood(model(test_x))\n",
    "    predictive_mean = trained_pred_dist.mean\n",
    "    lower, upper = trained_pred_dist.confidence_region()\n",
    "    \n",
    "f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "# Plot training data as black stars\n",
    "ax.plot(train_x, train_y, 'k*')\n",
    "# Plot predictive means as blue line\n",
    "ax.plot(test_x, predictive_mean, 'b')\n",
    "# Shade between the lower and upper confidence bounds\n",
    "ax.fill_between(test_x, lower, upper, alpha=0.5)\n",
    "ax.set_ylim([-3, 3])\n",
    "ax.legend(['Observed Data', 'Mean', 'Confidence'], bbox_to_anchor=(1.6,1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now our model seems to fit well on the data. It is not always possible to visually evaluate the model in high dimensional cases. Thus, now we evaluate the models with help of probabilistic metrics. We have saved predictive distributions from untrained and trained models as `untrained_pred_dist` and `trained_pred_dist` respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Negative Log Predictive Density (NLPD)\n",
    "\n",
    "Negative Log Predictive Density (NLPD) is the most standard probabilistic metric for evaluating GP models. In simple terms, it is negative log likelihood of the test data given the predictive distribution. It can be computed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Untrained model NLPD: 0.88, \n",
      "Trained model NLPD: -0.21\n"
     ]
    }
   ],
   "source": [
    "init_nlpd = gpytorch.metrics.negative_log_predictive_density(untrained_pred_dist, test_y)\n",
    "final_nlpd = gpytorch.metrics.negative_log_predictive_density(trained_pred_dist, test_y)\n",
    "\n",
    "print(f'Untrained model NLPD: {init_nlpd:.2f}, \\nTrained model NLPD: {final_nlpd:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean Standardized Log Loss (MSLL)\n",
    "\n",
    "This metric computes average negative log likelihood of all test points w.r.t their univariate predicitve densities. For more details, see \"page No. 23, Gaussian Processes for Machine Learning, Carl Edward Rasmussen and Christopher K. I. Williams, The MIT Press, 2006. ISBN 0-262-18253-X\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Untrained model MSLL: 0.82, \n",
      "Trained model MSLL: -0.49\n"
     ]
    }
   ],
   "source": [
    "init_msll = gpytorch.metrics.mean_standardized_log_loss(untrained_pred_dist, test_y)\n",
    "final_msll = gpytorch.metrics.mean_standardized_log_loss(trained_pred_dist, test_y)\n",
    "\n",
    "print(f'Untrained model MSLL: {init_msll:.2f}, \\nTrained model MSLL: {final_msll:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to calculate the quantile coverage error with `gpytorch.metrics.quantile_coverage_error` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantile 95% Coverage Error: 0.01\n"
     ]
    }
   ],
   "source": [
    "quantile = 95\n",
    "qce = gpytorch.metrics.quantile_coverage_error(trained_pred_dist, test_y, quantile=quantile)\n",
    "print(f'Quantile {quantile}% Coverage Error: {qce:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean Squared Error (MSE)\n",
    "\n",
    "Mean Squared Error (MSE) is the mean of the squared difference between the test observations and the predictive mean. It is a well-known conventional metric for evaluating regression models. However, it can not take uncertainty into account unlike NLPD, MLSS and ACE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Untrained model MSE: 0.20, \n",
      "Trained model MSE: 0.04\n"
     ]
    }
   ],
   "source": [
    "init_mse = gpytorch.metrics.mean_squared_error(untrained_pred_dist, test_y, squared=True)\n",
    "final_mse = gpytorch.metrics.mean_squared_error(trained_pred_dist, test_y, squared=True)\n",
    "\n",
    "print(f'Untrained model MSE: {init_mse:.2f}, \\nTrained model MSE: {final_mse:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean Absolute Error (MAE)\n",
    "\n",
    "Mean Absolute Error (MAE) is the mean of the absolute difference between the test observations and the predictive mean. It is less sensitive to the outliers than MSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Untrained model MAE: 0.38, \n",
      "Trained model MSE: 0.16\n"
     ]
    }
   ],
   "source": [
    "init_mae = gpytorch.metrics.mean_absolute_error(untrained_pred_dist, test_y)\n",
    "final_mae = gpytorch.metrics.mean_absolute_error(trained_pred_dist, test_y)\n",
    "\n",
    "print(f'Untrained model MAE: {init_mae:.2f}, \\nTrained model MAE: {final_mae:.2f}')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "d4d1e4263499bec80672ea0156c357c1ee493ec2b1c70f0acce89fc37c4a6abe"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
